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INTRODUCTION 
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This  report,  submitted  in  fulfillment  of  Data  Item  A001 , Contract  No. 
DAAG53-76-C-0127,  presents  the  results  of  a study  performed  by  DBA  Systems,  Inc. 
investigating  the  applicability  of  array  algebra  digital  terrain  modeling  and  data 
compaction . For  the  purpose  of  this  study  the  data  were  assumed  to  be  collected  along 
epipolar  lines. 

During  the  course  of  this  study  two  options  for  converting  the  collected  data 
into  regularly  spaced  terrain  elevations  were  evaluated.  The  first  of  these  is  to  perform 
an  "Array  Prediction",  or  functional  representation,  of  the  data  directly  in  the  epipolar 
coordinate  frame  (Section  3).  This  approach  allows  for  data  compaction  and  subsequent 
evaluation  at  uniform  intervals  in  a gridded  "map"  coordinate  system.  The  second  option 
considered  during  this  study  applies  the  principle  of  array  algebra  in  a piecewise  "trans- 
location algorithm"  (Section  4).  In  this  method  the  non-gridded  epipolar  coordinates 
are  first  converted  (translocated)  to  regularly  spaced  elevation  data  and  then  subsequently 
compacted  using  the  methods  of  "array  prediction". 

During  this  study  it  was  found  that  efficient  direct  array  prediction  algorithms 
could  be  derived  and  that  the  explicit  formation  of  the  regression  coefficients  could  be 
avoided.  Using  these  algorithms  the  array  requirement  for  rectangularly  spaced  data 
was  removed;  thereby  allowing  the  boundaries  of  the  area  to  assume  arbitrary  shapes. 

In  delineating  these  methods  several  examples  of  direct  prediction  coefficients,  or 
"k-vectors",  are  derived.  In  typical  cases,  the  knowledge  of  ten  (10)  such  coefficients 
replaces  the  data  handling  as  conventionally  required  for  the  explicit  formulation  and 
solutions  of  normal  equations. 

The  algorithms  developed  or  translocation  demonstrate  that  between  2 and  10 
scalar  multiplications  are  required  for  each  point  manipulated.  During  the  course  of 
this  effort  the  various  alternatives  for  performing  the  operation  were  examined. 
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In  addition  to  analyzing  the  mathematical  equations  required  for  terrain 
compaction,  the  computational  requirements  were  analyzed  for  both  sequential 
and  parallel  processors.  Based  upon  this  analysis  it  ' ppears  that  the  associative 
array  processor  is  ideally  suited  for  processing  the  derived  algorithms.  This  is  due 
to  the  feature  of  "sweeping  while  scanning",  i.e.  processing  the  algorithms  first  in 
one  direction  and  immediately  using  the  results  in  the  other  coordinate  direction. 

It  appears  that  only  minimal  storage  is  required  to  generate  the  one  sweep  characteristic 
for  the  algorithms,  thereby  avoiding  multiple  read-in  (fetch)  or  read-out  (store)  cycles. 

The  remainder  of  this  report  describes,  in  detail  the  results  of  this  effort. 
Section  2 provides  some  relevant  background  on  the  principle  of  array  algebra  and 
illustrates  its  application  using  simple  functional  examples.  Sections  3 and  4 put 
forth  the  detailed  mathematical  developments  and  analysis  of  the  algorithms  for  direct 
prediction  and  translocation.  The  final  section  (5)  presents  the  conclusions  and 
recommendations  based  upon  the  study. 
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BASIC  PRINCIPLES  OF  ARRAY  ALGEBRA 
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This  section  is  devoted  toon  exposition  of  the  basic  principles  of  Array  A Igebra , 
centered  around  the  notions  of  array  multiplication  and  array  least  squares  regression. 

The  first  subsection  introduces  Array  Algebra  in  the  simplest  situation  of  two  dimensions, 
and  contains  a numerical  example  to  illustrate  the  ideas.  The  second  subsection  presents 
the  theory  for  higher  dimensions.  The  third  subsection  discusses  the  savings  in  computations 
that  array  algebra  solutions  to  linear  equations  produce  compared  to  conventional  least 
squares  techniques. 

2 . 1 Array  Algebra  in  Two  Dimensions  with  a Numerical  Example 

Let  m and  n be  positive  integers.  By  a 2-dimensional  array  of  size  (n  x m)  is 
meant  a collection  of  real  (or  complex)  numbers,  indexed  by  two  indices  t and  j,  which 
assume  integer  values  between  1 and  n and  m,  respectively.  An  individual  number  of 
the  collection,  called  an  element  of  the  array,  is  denoted  by  the  form  xn  , and  the 
array  itself  is  abbreviated  [x, , ] . 

As  an  ordinary  matrix  is  also  a collection  of  numbers  indexed  by  two  indices, 
the  separate  notion  of  2-dimensional  array  may  appear  superfluous.  It  will  be  seen  below, 
however,  that  the  roles  of  arrays  and  matrices  are  different. 

MtntceA  witl  be  denoted  by  capital*  fin om  the  beginning  o<(  the 
alphabet,  aAAay*  by  capitals  fan  om  the  end. 

Although  a 2-dimensional  array  has  the  same  structure  as  a matrix,  its  role  in 
mathematics  is  in  several  ways  like  that  of  a vector.  In  particular,  two  arrays  X [ xf  s J 
and  Y [_y,  j ] of  the  same  size  (n  x m)  may  be  added  together  to  form  a new  (n  x m)  array 
as  follows: 

X + Y = [xfJ  + yn]  . 
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This  operation  is  called  array  addition . Also,  if  c is  a scalar,  that  is,  a real  (or 
complex)  number,  then  c and  the  array  X may  be  multiplied  to  form  a new  (nxm)  array 
as  follows: 

cX  = [ex,  3 .1  . 

This  operation  is  called  array  scalar  multiplication. 

The  two  array  operations  defined  above  satisfy  the  same  rules  as  do  the 
corresponding  operations  on  vectors.  Because  of  this,  the  set  of  all  arrays  of  a given 
size,  together  with  the  two  array  operations,  have  the  structure  of  a vector  space. 

For  given  size  (nxm),  this  space  will  be  denoted  , and  will  be  called  (nxm) 
array  space . 

There  is  a third  type  of  operation  on  an  array  space  Rn,n  which  resembles 
the  operation  of  multiplying  a vector  by  a matrix  to  produce  a new  vector. 

Let  X = [x,  j ] be  an  (n  x m)  array;  let  A = [ah  k ] be  a (pxq)  matrix . If 
q=n,  an  operation  denoted  oa  on  A and  X that  yields  a new  (pxm)  array  is  defined 
as  follows 1 : 

A °i  X = J = Cyhl  ]• 

[p,n]  (n,m)  (p,m) 

The  new  array  is  indexed  by  h and  j,  in  that  order  and  equates  to  the  hthrow  of  A times 
the  3 th  column  of  X . 

If  q = m,  a second  operation  denoted  o2  on  A and  X that  yields  a new  (nxp) 
array  is  defined  as  follows'  : 

ID 

A o3  X = L £ ahlf  x1k  ) Lylh]  ■ 

If  = 1 

(p,m)  (n,m)  (n,p) 

The  new  array  is  indexed  by  i and  h , in  that  order  and  equates  to  the  hth  row  of  A times 
the  i th  row  of  X . 

1 The  operation  K ox  li,  denoted  by  a MipesiAcAtpt  A1  in  Vauhala  ( 7 974 ) 

2 Denoted  A 2X  in  Rauh&ia  ( 1974 ) 
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Both  of  the  above  operations  will  be  called  array  multiplication.  They  havo 
the  property  of  being  linear  operations.  That  is, 

A Oj  (cX  + dY)  = cA  Oj  X + dA  o;  Y, 


A o2  (cX  + dY)  = cA  Oa  X + dA  Og  Y. 

These  equalities  follow  directly  from  the  definition  of  the  operations,  in  the  same  way 
as  do  the  corresponding  operations  on  vectors. 

Three  further  important  equalities  are  as  follows: 


(AB)  o1  X = A o,  (B  o1  X), 
(AB)  o2  X = A o2  (B  o3  X), 


where  (AB)  = matrix  product  of  A and  B 


° A 0l  (Bo2  X)  = Bo,  (A  Oj  X). 

These  equations  are  derived  by  observing  that  if  A = [ahk  J,  B=  [b^  ] and  X = [xl3]  , 
then  both  sides  of  the  first  equation  may  be  written  as 

S ahIc  by.,  x,  i J , 

l k 

and  both  sides  of  the  second  equation  may  be  written  as 

^ ahk  bk  3 Xjj  ] . 
lc  3 

For  the  third  equation,  observe  that 

^ ah  1 ^ J X1  J ^ J ^ °M  X1  J • 

<3  3 ! 


A fourth  operation  on  arrays  is  the  dot  product,  which  is  similar  to  the  dot  product 
for  vectors.  Namely,  if  X = [xfj  J and  Y=  [y,j]  are  both  arrays  of  size  (nxm),  then 

X • Y = L L xl3  yM  ; 

1 3 « 

and  is  a number  equal  to  the  sum  of  the  products  of  all  corresponding  elements  in  the  two 
arrays  1 . 

1 Ve.no ted  RS ( X* V ) in  Rauhala  [1974] 
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The  value  of  a dot  product  is  a scalar.  The  dot  product  operation  is  only  defined  for 
arrcys  of  the  same  size.  If  complex  numbers  are  involved,  the  above  definition  must  be 
changed  to 

^ = ^ xu  yu  ' 

where  the  bar  denotes  complex  conjugation. 

Recall  that  if  A = [at]k  ] is  an  (nxm)  matrix,  then  the  transpose  A1  of  A is  the 
(mxn)  matrix  whose  (k,h)th  entry  a‘k  h is  ahk  . The  transpose  operation  is  related  to  the 
dot  product  for  vectors  as  follows.  If  v is  an  m-vector,  and  u is  an  n-vector,  then 

(Av)  • u = v • (AT  u) . 

There  are  corresponding  equations  for  the  array  dot  product  and  the  array 
multiplication  of  matrices  and  arrays:  if  A is  an  (nxm)  matrix,  and  X and  Y are  (mxp) 
and  (n  x p)  arrays  respectively,  then 

(A  Oj  X)  • Y = X-  (Ar  0l  Y). 

This  equality  follows  from  the  fact  that 

2 2 £ an  xix  )yIk  = 22  Xjk  £ an  y,k  )• 

1 k J J k l 

If  now  X and  Y are  (pxm)  and  (pxn)  arrays  respectively,  then 
(A  o3  X)  • Y = X • (AT  os  Y). 

This  follows  from 

SL  (I  ak  , xn)y,k  = £ £ x j (2  ak } y,k ) . 

1 V .1  1 .1  tr 

In  case  complex  numbers  are  involved,  the  transpose  operation  must  be  replaced 
by  the  complex  transpose*;  namely,  A*  is  the  matrix  whose  ( J, s)th  element  a!(  equals 
rhe  complex  conjugate  a, , of  the  (l,t)th  element  a, s of  A . 
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The  numerical  example  of  this  subsection  is  a 2-dimensional  array  least  squares 
regression . 

Let  Pj  and  p2  each  be  a finite  set  of  parallel  lines  in  the  Euclidean  plane  E2 , 
but  such  that  the  lines  of  px  are  not  parallel  to  the  lines  of  pa  . The  set  of  all  points 
that  are  intersections  of  lines  pT  and  p2  is  called  a bounded  grid  in  E2  . 

In  this  example,  let  p2  consist  of  the  lines  x = 0,  x = 1 , x = 4,  x = 5;  p2 

consist  of  the  lines  y = -4,  y - 1 , y = 2 . The  associated  bounded  grid  G therefore 
consists  of  the  12  pairs  of  points  (x,y)  such  that  x equals  either  0,1,4  or  5,  and  y equals 
either  -4,  1 or  2 . 

Associated  with  G is  a pair  (i  , j)  of  indices,  1 S i S 4,  1 S j S 3,  such  that  the 

(l  ,l)th  element  of  G is  the  intersection  of  the  ith  line  of  px  and  the  j th  line  of  p2;  the 

set  of  parallel  lines  is  itself  ordered  either  from  left  to  right,  or,  if  horizonal,  from 
bottom  to  top. 

To  the  (i,  t)th  point  of  G,  let  there  correspond  a number  u(J  such  that  the  (4x3) 
array  U = [u<f]  is  as  follows 

2 -1 

1 1 

0 2 

-1  0 

The  elements  uf , of  U can  be  considered  as  values  of  an  empirically  given  function  at  the 
grid  points  (x,  ,yj)  of  correspoinding  index  pairs  (l , j)  • Analogously,  if  the  pairs  (x,y) 
of  the  grid  are  latitudes  and  longitudes  of  points  on  the  Earth,  the  u] , could  represent 
radial  elevations  above  a regular  figure  such  as  a sphere  or  spheroid  fitted  to  the  surface 
of  the  Earth . 


0 

1 

-1 

-3 


-7- 


The  problem  is  to  determine  the  coefficients  wh  in  the  math  model 

z = s s wt.ic  fhW  9 if  (y) 

h k 

which  "best"  fit  the  empirical  values  U . Two  finite  sets  of  functions  f: ,fp  and 
gT  , . . .,g5  of  a single  real  variable  are  given,  such  that  the  domains  of  the  fh  include 
all  of  the  first  coordinates  of  the  grid  points  in  G,  and  the  domains  of  the  gk  include 
all  of  the  second  coordinates  of  the  grid  points  of  G.  Then  real  values  for  variables 
wh]t  , lShSr,  ISk-s  are  to  be  determined,  so  that  the  sum 

i £ (£  £ fh(x,)9,(y3)whk  -un)I * 3 

1 = 1 J — 1 h=  1 k=  1 

is  a mi nimum  . 


In  the  present  example , let  r = 3 and  s = 2 . A Iso , assume  fh  (x)  = xh  _1,  and 
gk(y)  = yk_1  • The  functions  fh  and  gk  and  the  grid  G give  rise  to  two  matrices 
A = [fh (x, )]  and  B = [gk  (y_, )1  of  sizes  (4x3)  and  (3x2)  respectively.  Namely 


1 0 0 

1 1 1 

1 4 16 

1 5 25 


B 


1 

1 

1 


Form  the  (3x2)  array  W = [whk  ] . Then  using  the  array  dot  product,  the  above 
sum  of  squares  may  be  expressed  as 


(B  os  A 0l  W - U)  • (B  o8  A o,  W - U). 

Tfie  objective  Is  to  ex.ps.eM  the.  solution  value  fax  the  axxaij  U In  this 
minimization  a*  an  axxay  pxoduct  Involving  the  axxay  U.  This  Is  called  axxay 
least  squaxes  Xegnesslon . 

I -i 

For  any  (nxm)  matrix  C of  rank  m,  for  which  mSn,  set  C = (CTC)  CT  , 
or,  if  C contains  complex  numbers,  C^=  (C*C)  C*.  Then  as  shown  below,  the  solution 
value  for  W is 

W0  = B^  o2  A^  Oj  U . 
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Recall  that  two  vectors  x and  y in  a vector  space  are  orthogonal  if  their  dot 
product  x *y  equals  0.  Using  the  analogous  notion  of  orthogonality  in  the  array  spaces,  W0 
may  be  characterized  by  the  condition  that  B os  A o,  W0  - U is  orthogonal  to  B oE  A Oj  W 
for  all  arrays  W in  Rr  s = R3  2 . Geometrically  this  means  that  the  array  B o?  A Oj  Wn  - U 
is  orthogonal  to  the  subspace  of  all  B o2  A ox  W , and  so  that  B oE  A o(  W0  is  the  element 
of  this  subspace  closest  to  U . It  must  therefore  be  shown  that 

(B  Oj  A Oj  W)  • (A  Oj  B o,  B Og  A ox  U - U)  = 0 

for  all  W in  R3  2 . This  equation  may  be  rewritten  as 

W • (BT  o2  AT  Oj  A Oj  B Oj  B^  o2  A^  ox  U - BT  o2  AT  ox  U)  = 0. 

But  ^ ^ 

BT  o2  AT  ox  A Oj  B Oj  B o2  A ox  U 

= (BTBB^)  o2  (ATA  A")  ox  U 

= BT  o2  AT  ox  U . 

Thus  the  right-hand  side  of  the  above  dot  product  equals  the  zero  array,  so  the  dot 

-t  t 

product  itself  equals  0 for  all  W in  R3  2 . Hence  W0  = B o2  Aq  ox  U . In  order  to 

' i l 

compute  the  value  W0  of  W for  the  example,  it  is  therefore  necessary  for  form  A ,B  , 

'L  'L  £ 

and  to  evaluate  B o2  A ox  U . The  (3x4)  matrix  A is  as  follows: 


' .86765 

.22059 

-.22059 

.13235 

-.77206 

.53676 

.71324 

-.47794 

.12500 

-.12500 

-.12500 

.12500 
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The  (2x3)  matrix  B equals 

.27417  .35484  .37097 

-.17742  .06452  .11290  ‘ 

f 'L  £ £ 

Using  these  values  for  A and  B , the  array  W0  = B o2  A o:  U is  found 

to  be 

.25546  -.01305  ' 

1.17161  .05159  • 

-.28831  .01008 


2 .2  Array  Algebra  in  Higher  Dimensions 

To  simplify  the  writing  of  arrays  in  higher  dimensions,  the  following  notion 
of  index  vector  is  introduced.  Let  i,  ...jn  be  a sequence  of  integer  valued  indices 
such  that  i takes  values  1 to  nk  . Form  the  index  vector  J = (jt,...,1n).  Then  J 
takes  as  values  the  vectors  with  integer  coordinates  in  the  hyperrectangle  Q which 
consists  of  all  vectors  (xj ,..  .,xN)  such  that  l-xk  5nk  , iSkSN. 
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An  N -dimensional  array  of  size  (na  x n2  x . . .x  nN)  is  a collection  of 
real  (or  complex)  numbers,  indexed  by  N indices  or  more  compactly, 

by  the  index  vector  J - (^,...,3  ,)•  The  index  vector  J takes  as  valueS  vectors 
with  integer  coordinates  in  a hyperrectangle  Q . An  array  element  is  written  in 

the  form  Xj , and  the  array  itself  is  abbreviated  [Xj  1 . 

As  in  the  case  of  two-dimensions,  the  set  Rn  , of  all  arrays  of  size 

(nxx...x  nN)  form  a vector  space,  called  an  array  space.  The  two  operations  of  array 
addition  and  scalar  multiplication  are  defined  by 

[xj  ] + Lyj  ] = [xj  + Yj  3 / 

and 

c[Xj]  = [cXj]  . 

A matrix  A = [atl  ] and  an  N-dimensional  array  [xj  ] can  be  multiplied 
together  in  N different  ways,  each  called  an  array  multiplication1.  Namely,  if  A 
has  size  (mxn)  and  [xx  1 has  size  (nx  x...xnj,  then  if  n=  nk  let 


A ok  X = (2  a,  j Xj) . 
l.  k 


AokX  has  size  (nx  x. . .x nk _x  xmxnk+1x, 
/ 3k+  1 ' ' • • '5n 


.x nN ),  and  index  vector 


As  for  the  case  of  two  dimensions,  array  multiplication  in  N dimensions  has  the 
following  properties: 


A ok  (cX  + dY)  = cA  ok  X i dA  ok  Y, 

(AB)  ok  X = A ok  B ok  X, 

A ok  B oh  X = B oh  A ok  X , 

for,  yf  h.  Fork  = h , the  second  relation  shows  that  the  third  relation  holds  only  for 
A and  B such  that  AB  = BA . 


1 Tunmed  R matAlx  muJUlpllcdtAoni  in  Rauhala  (1974) 

2 l Idnotzd.  ai  A*  X in  Rauhaia  ( 1974) 
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The  dot  product  of  two  arrays  X and  Y of  the  same  size  is  defined  by 

X • Y = 2 Xj  yj  , 

j 

where  j varies  over  the  vectors  with  integer  coordinates  in  the  hyperrectangle  Q. 

To  formulate  the  problem  of  N-dimensional  array  least  squares  regression,  the 
concept  of  a bounded  grid  in  N-dimensional  Euclidean  space  EN  must  first  be  defined. 

Let  nk  be  a positive  integer  for  15*  SN.  For  each*  , let 

be  real  numbers . The  set  G of  all  vectors  in  EN  of  the  form 

*•  n * , . 

X = (x(l) x,  )T  , where  15  jk  5mk  , 15*  5 N,  is  called  a rectangular  bounded 

J J1  3 M 

grid  in  E*,  . 

Let  B be  a nonsingular  matrix  of  real  numbers  of  size  (NxN).  Then  the  set 
D = BG  of  all  vectors  Bxj  , where  Xj  is  a vector  in  the  rectangular  bounded  grid  G, 
is  called  a bounded  grid  in  EN  . 

Let  e; , . . . ,eN  be  the  standard  orthonormal  basis  for  EN  , so  e^  (1 ,0, ...  ,0), , 

eN  = (0, . . .,0,1).  Then  Bet  ,...,  BeN  Is  also  a basis  in  EN , but  it  is  in  general  not 

rectangular . Set  bk  = B ek  . If  y - ^ , • • • r Yu ) = 2 yk  e*  , then  By  - B(L  yv  ek  ) 

£ yk  Bek  = £ yk  bk  . Therefore  the  elements  of  the  bounded  grid  D may  be  expressed 

V Mu 

as  2j  x j bk  . 

Ic  J k 

Let  U = [ujl  be  an  N-dimensional  array  of  size  (nx  x. . .xnj  whose  elements 
are  in  practice  values  of  an  empirically  defined  function  at  the  points  Bj  of  the  bounded 
grid  D in  EN  . 

For  each  * , 15*  5N,  let  m„  be  a positive  integer,  and  let  fk  , , 1-  -mk  / 

K k 

be  functions  of  a single  real  variable,  such  that  for  each  index  vector  j , the  v th 
coordinate  coefficient  x^  of  Bxj  with  respect  to  the  basis  b1 , . . . ,bN  belongs  to  the 

3 * (k)  _ 

domains  of  fk  , , . . . ,fk  „ ^ , 1 5 * 5 N . Let  A k be  the  (nk  x mk  ) matrix  rfMk(x,k)j  • 
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The  problem  of  array  least  squares  regression  is  to  find  the  array  W of  size 
(m;  x . . .v  mN ) such  that  the  dot  product 

(Aj  ox  A2  o2  • • • AN  o^  W — U)  (A:  o^  A2  o2  » . . An  on  W — U) 

is  a minimum . 

By  an  argument,  based  on  orthogonality,  entirely  similar  to  the  case  of  two 
dimensions,  the  solution  value  of  W is 

W0  = aJo,  A2  o2...  A^  on  U, 

where 

aJ  - (a;  a,)-1  a;. 

In  case  the  elements  of  U are  measured  quantities,  and  the  variances  of  all 
the  measurements  are  the  same  and  equal  to  the  constant  crs  , then  the  maximum  likelihood 
estimate  of  U of  the  form  A1  o:  A2  o2  ...  AN  oN  W is  the  value  of  this  expression  for 
which  W equals  the  W0  above.  This  follows  readily  from  the  fact  that  the  covariance 

N 

matrix  of  the  observations  is  a3  I , where  I is  the  identity  matrix  of  order  7 T nk  . 

k=  l 

2.3  Computational  Advantages  of  Array  Algebra  and  Array  Regression,  and 
Suitability  for  Programming  on  a Parallel  Array  Processor 

The  role  of  array  algebra  in  computational  applications  consists  of  being  a 
formalism  for  an  extremely  efficient  method  of  solving  systems  of  multilinear  equations 
in  which  certain  of  the  sets  of  quantities  involved  are  associated  with  a bounded  grid  of 
points  in  some  higher  dimensional  Euclidean  space.  Both  the  method  of  solution  and  the 
formalism  of  array  algebra  arise  out  of  this  association. 


By  the  use  of  the  operation  of  array  multiplication,  the  amount  of  computation 
required  for  the  solution  of  such  a multilinear  system,  e .g . , in  the  form  A1  oT  . . .AN  oN  U , 
where  all  the  matrices  Ax , . . . ,AN  have  the  same  size  (nxn),  is,  when  done  by  sequential 
processing,  proportional  to  the  first  power  of  the  number  nN  of  parameters,  and  not  to  the 
third  power(nN)3as  in  the  conventional  linear  solution. 

The  array  solution  method  may  also  be  compared  with  another  fast  algorithm 
for  solving  systems  of  linear  equations  of  a special  type,  namely  the  Fast  Fourier  Transform. 
For  this  algorithm,  the  amount  of  computation  required  is  proportional  to  the  product  of 
the  number  of  parameters  times  the  logarithm  of  the  number  of  parameters.  Moreover, 
this  level  of  efficiency  is  only  attainable  when  the  number  of  parameters  is  a power  of  2. 

If  A k is  an  (mk  x nk  ) matrix,  and  X is  an  nxx . . .x  nN  array,  then  on  a parallel 

array  processor,  the  operation  A ok  X requires  mk  nk  parallel  multiplications.  Thus  array 

# 1 

products  of  the  form  A1  ox  . . . AN  oN  X or  Ax  Oj  • • • AN  oN  U , in  which  the  array  U has 
size  mx  x. . .xm„  , require  a sequence  of  N subsequences  of  m^  , . . .,  mNnN  parallel 
multiplications  respectively.  This  shows  that  array  algebra  operations  and  array 
regression  are  eminently  suited  for  implementation  on  a parallel  array  processor. 
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3. 


APPLICATION  OF  ARRAY  LEAST  SQUARES  REGRESSION  TO  DIGITAL 
TERRAIN  DATA  COMPACTION  IN  THE  EPIPOLAR  COORDINATE  SYSTEM 


The  basis  for  the  application  of  array  least  squares  regression  to  the  compaction 
of  digital  terrain  data  in  the  epipolar  coordinate  system  is  that  the  elevation  observations 
are  made  at  a set  of  points  which  in  the  epipolar  coordinate  plane  form  a grid.  This 
implies  that,  in  the  terminology  of  Section  2,  the  observations  form  a 2-dimensional 
array.  The  functions  that  are  used  in  the  regression  have  a special  form  which  is 
appropriate  for  terrain  modelling;  these  functions  are  described  in  subsection  3.1.  The 
structure  of  the  compaction,  the  coefficients  for  the  regression,  and  the  array  algebra 
formula  that  expresses  these  coefficients  in  terms  of  the  observations  are  discussed  in 
subsection  3.2.  Subsection  3.3  addresses  the  problem  of  specifying  the  values  of 
certain  parameters  that  determine  the  individual  functions  of  the  regression.  These 
parameters  are  not  solved  for,  but  are  specified  in  advance  on  the  basis  of  the  general 
structure  of  the  terrain  being  modelled  and  the  density  of  the  elevation  observations 
made  by  the  epipolar  scanner.  The  recovery  from  the  compaction  of  the  elevations  at 
all  of  the  observation  points,  for  input  to  the  contour  plotter,  is  described  in  sub- 
section 3.4. 

3 . 1 Functions  Occurring  in  the  Array  Least  Squares  Regression  for  Digital  Terrain 

Data  Compaction 

Let  U = [utJ]  be  the  2-dimensional  array  of  elevations  observed  by  an  epipolar 
scanner  at  a set  of  observation  points  which  in  the  epipolar  coordinate  system  form  a 
bounded  grid  D.  By  Section  2.2,  there  is  a nonsingular  linear  transform  B and  a 
rectangular  bounded  grid  G such  that  D = BG . If  bt  = B(1,0)T  and  b?  = B (0 , 1 )'  , then 
for  each  (x,y)  in  the  bounded  grid  G,  xbj  + yb?  is  an  element  of  the  bounded  grid  D. 

Let  x;  , . . . ,xN  and  yx  , . . . ,y„  be  real  numbers  such  that  G consists  of  the  vectors 
(*,  » y . );  then  D consists  of  all  x,  bx  + y , b2  , for  1 5 i = N and  1 § i 5 M . 
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According  to  Section  2,  the  problem  of  array  least  squares  regression  in  two 
dimensions,  as  applied  to  the  array  U of  observations  at  grid  points  in  D,  is  to  determine 
the  coefficients  a,,  , in  a linear  combination 

2 2 ahk  fh(x)  gk  (y) 

h k 

such  that  the  sum  of  squares 

2 2 (£  2 ahk  fh(x, ) 9k  (/3)  - Un)8 

1 J h k 

is  a minimum . 

The  form  of  the  functions  fh  and  gk  to  be  considered  for  the  present  application 
to  digital  terrain  data  compaction  is  as  follows.  Let  c be  a real  number,  c<  0.  Then  the 
function  exp  (cx2)  is  called  a Gaussian  error  function,  and  the  functions  fh  and  gk  will 
be  translates  of  functions  of  this  form.  That  is,  there  are  constants  ch  and  sh  , dk 
and  tk  such  that 


and 


fh(x)  = exp  (c„  (x  - shf), 


9k  (y)  = exp  (dk  (y  - tk  J3). 


If  IS  uSH  and  1 Sk  S K,  then  the  points  (sh  ,tk  ) are  called  node  points  for  the  regression. 


1 


Gaussian  error  functions  take  the  value  1 at  the  origin  0,  and  are  symmetric 
about  the  origin  because  the  argument  in  the  exponential  is  squared.  They  also  decrease 
rapidly  to  0 as  the  argument  grows  in  magnitude,  the  rapidity  being  greater  the  more 
negative  the  constant  c. 


The  motivation  for  the  choice  of  translates  of  Gaussian  error  functions  for  the 


f 


f 
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regression  is  that  for  the  prediction  or  estimation  of  the  elevation  at  a point  with 
coordinates  Xq  bx  + y0  b3 , by  the  evaluation  of  the  functional  model 

F(x,y)  = 2 2 ahk  fh(x)  gk  (y) 

h k 

at  ,y0 ),  the  values  fh  ) and  gk  (y0)  for  node  points  (sh,tk  ) not  very  near  to  (x,  ,y0  ) 
are  negligibly  small,  so  that  F(x0  ,y0  ) is  adequately  expressed  by  a sum 


2 2 ahk  fh  (Xq  ) 9k  (y0 )/ 

h k 

where  h and  k take  only  a small  range  of  values  corresponding  to  node  points  (st  ,tk  ) 
very  close  to  (>v,  ,y0). 


A more  detailed  discussion  of  the  origin  and  suitability  of  the  above  choice  of 
translates  of  Gaussian  error  functions  for  the  regression  is  contained  in  Kraus  and  Mikhail 
[2  ]. 


3.2  The  Nature  of  the  Compaction  and  the  Array  Algebra  Formula  for  the 
Coefficients 

The  functional  model 


F(x,y)  = 2 2 ahk  fh  (x)  gk  (y) 

h k 

being  the  one  used  to  estimate  terrain  elevations,  it  follows  that  the  compaction  of  the 
terrain  data  consists  of  the  following: 

(a)  The  set  of  node  points  (sh  ,tk  ),  and  the  coefficients  ck  and  dk  , !5n5H, 

1 5 y $ K . These  are  selected  on  a basis  to  be  discussed  below . 

(b)  The  set  of  coefficients  ahk  , lStiSH,  l^kSK.  These  are  computed  using 
the  expression  given  in  Section  2 for  the  solution  to  the  2-dimensional  array  least  squares 
regression  problem. 


Relative  to  the  computation  of  the  coefficients  ahk  , it  will  happen  in  practice 
that  the  number  NM  of  observation  points  is  very  large,  as  may  be  the  number 
HK  of  coefficients  ahk  . Nevertheless  in  the  array  algebra  determination  of  the 
coefficients  ahlc  , the  only  matrices  that  need  to  be  inverted  are  the(HxH)matrix  A*AX 
and  the(KxK)matrix  A2A2  , where 

A1  = [fh(x,)],  of  size  (NxH), 

and 

A2  = [9*  (ys )],  of  size  (MxK). 

; _1  t _1 

Setting  At  = (A^AX)  Ax  and  A2  = (A^Ag)  A\  , then  in  the  notation  of  Section  2, 
the  array  [ahk  "|  is  equal  to 

Ai  °i  *£  U . 

This  formula  is,  by  the  discussion  at  the  end  of  Section  2.2,  valid  for  measured 
quantities  for  which  all  the  associated  variances  are  taken  to  have  the  same  values. 

A slightly  more  general  situation  may  be  accommodated  by  the  array  algebra; 
this  is  the  case  in  which  the  covariance  matrix  for  all  the  observations  is  expressible  as 
a tensor  product  Cx  £>C2  . Here  Cx  = [c,^  ] is  of  size  (NxN),  C2  = [cj'j  ] is  of  size 
(MxM),  and  Cx  ® Cs  is  by  definition  equal  to  the  (NMx  NM)  matrix  [Cx  c("3  ] . In 
this  case  the  array  of  coefficients  [ah)c  ] equals 


((A;  C"1A1)"1AI  C^1)  ox  ((A l C2_1A2)_1A;  C'2l)  os  U . 


In  practice,  the  matrices  Cx  and  C3  will  be  diagonal,  so  their  inversions  are  performed 
by  inverting  their  scalar  diagonal  elements  term-by-term. 


e 

t 
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Also  in  practice,  an  empirically  determined  covariance  matrix  may  be 
approximated  more  or  less  well  by  a positive  definite  tensor  product  matrix.  This 
could  be  used  in  the  array  formula  above  for  the  regression  coefficients  to  produce 
values  for  these  which  are  better  in  the  sense  that  their  variances  are  smaller  than 
would  otherwise  h*  f the  covariance  matrix  of  the  observations  was  tacitly  assumed 
to  be  a positive  scalar  multiple  of  the  identity  matrix,  the  scalar  being  realistically 
large. 

3 .3  The  Specification  of  the  Parameters  in  the  Gaussian  Error  Functions 

The  functional  model  of  the  terrain  elevations  is  a linear  cotnbination  of 
products  of  pairs  of  functions  fh  and  gk  , each  such  function  being  selected  as  a 
translate  of  a Gaussian  error  function,  hence  of  the  mappings 

x -»exp  (cfc-s)2),  y -*  exp  (d(y  - t)3), 

where  c,d<0and  s,tare  parameters.  This  subsection  addresses  the  problem  of  specifying 
the  values  of  the  parameters  c and  d,  and  s and  t in  all  of  the  functions  f„  and  gk 
occurring  in  the  functional  model. 

Two  considerations  influence  the  choice  of  the  parameters.  These  are  computational 
feasibility  and  efficiency,  and  the  amount  of  undulation  of  the  terrain  being  modelled,  i.e., 
whether  it  is  level  to  broadly  rolling,  or  very  rolling  to  broken.  These  two  considerations 
are  interconnected,  and  will  be  discussed  simultaneously. 

To  avoid  ill-conditioning  in  the  matrices  to  be  inverted  in  the  determination  of 
the  coefficients  ahk  , it  is  necessary  that  for  a given  spacing  of  the  node  points  the 
magnitudes  of  the  parameters  ch  and  dk  in  fh  and  gk  respectively  are  not  too  large  or 
too  small . If  the  coordinate  system  is  so  chosen  that  the  mean  distance  between  adjacent 
node  points  is  1,  then  the  ch  and  dk  should  be  not  more  than  Tn(0.8)  = - .223,  and  not 
less  than  ln(0.2)  = -1 .609.  Within  these  limits,  all  the  ch  are  assigned  a single  value  c, 
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and  all  the  dk  are  assigned  a single,  possibly  different,  single  value  d.  These 
assignments  are  made  with  a view  to  simplifying  computations.  The  actual  values  of 
c and  d,  within  the  limits  set  above,  will  be  farther  from  0 for  less  rolling  terrain,  and 
closer  to  0 for  more  rolling  terrain. 

As  will  be  shown  in  3.4,  computational  advantages  occur  if  all  the  node  points 
are  also  observation  points.  Then  the  node  points  will  form  a subgrid  of  the  grid  of 
observation  points.  This  is  especially  true  if  the  node  point  grid  is  obtained  from  the 
observation  point  grid  by  deleting  every  other  row  and  every  other  column.  This  is  the 
situation  that  will  be  considered  here.  The  compaction  will  then  consist  of  the  at  most 
(N+1)(M+l)/4  coefficients  ahk,  the  values  of  the  parameters  ch  , dk  , sh , tk  being 
specifiable  by  only  a few  words. 

It  is  clear  that  because  the  spacings  of  the  observation  points  and  the  node 
points  are  closely  related,  and  because  the  spacing  of  the  node  points  and  the  values  of 
the  parameters  ch  and  dk  are  also  related,  the  observation  point  spacing  must  be  chosen 
in  a particular  way.  Namely,  the  spacing  of  the  observation  grid  points  is  such  that 
relative  to  the  resulting  subgrid  of  node  points,  values  of  ch  and  dk  can  be  chosen  in 
the  interval  from  -1 .609  to  - .223,  as  described  above,  so  that  undulations  of  the  terrain 
are  adequately  modelled  by  the  function 

F(x,y)  = £ £ ahk  fh  (x)  gk  (y). 

h k 

3 .4  Recovery  of  Observation  Point  Values  from  the  Compaction 

In  the  most  general  situation,  in  which  the  node  points  do  not  necessarily 
form  a subgrid  of  the  grid  of  observation  points,  elevation  values  at  grid  points  are 
recovered  from  the  compaction  in  the  form  of  the  estimates  obtained  by  evaluating 
the  functional  model  at  the  observation  grid  points.  If,  as  in  subsection  3.2, 
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Ai  = [f„  (x,  )]  and  As  = [gk  (yj)]  , where  (x,  ,yj),  1 S iS  N,  IS  are  the 

observation  grid  points,  then 

Ai  °1  A2  °s  C°hk  3 

is  the  array  of  observation  point  estimated  elevations. 

Because  of  the  choice  of  the  functions  fh  and  gk  , the  matrices  Aq  and  A3  may 
in  practice  be  replaced  by  banded  matrices  Aq  and  Ag  , obtained  from  Aq  and  Ag  by 
replacing  all  of  the  negligibly  small  elements  in  them  by  0.  In  a single  column,  these 
will  be  all  elements  sufficiently  far  from  the  diagonal;  i.e.,  corresponding  to  numbers 
x,  or  y3  which  are  far  enough  away  from  the  sh  or  tk  for  that  column  respectively,  so 
that  fh  (x, ) = exp  (ch  (x,  - sh)2)  or  gk  (y,)  = exp  (dk  (y3  - tk  )s)  are  negligibly  small. 

In  practice,  for  this  general  situation,  the  array  of  elevations  at  the  grid  of 
observation  points  is  estimated  by 

Ai  °i  As  °2  ^ahk  "I  • 

If  the  bandwidth  of  A^  is  2ba  and  that  of  Ag  is  2b3  , then  (2b1  + l)NK  + 
(2bg+l)NM  scalar  multiplications  are  required  in  evaluation  of  the  NM  elevation 
grid. 

When  the  node  points  form  a subgrid  of  the  grid  of  observation  points,  as 
described  in  subsection  3.3,  then  by  a transformation  of  the  functional  model,  the  array 
of  coefficients  [atlk  1 may  be  replaced  by  a new  coefficient  array  [F(sp,tq  )T  of  values 
of  the  functional  model  at  the  grid  of  node  points  (s  , tq  ),  1 3 p 5 H , 1 5<i  3 K . This 
does  not  achieve  any  further  compaction,  but  does  reduce  the  amount  of  computation 
necessary  to  recover  the  elevations  at  all  the  observation  grid  points. 

Let  Bi  = [fh  (sp)Jand  B2  = [gk  so  that  B:  and  B3  have  sizes  (H  xH)  and 

(Kx  K)  respectively,  and  are  moreover  nonsingular.  Then 

[F(sp,tq)]  = Bj  oq  Bg  Og  [ahk  1; 
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!-ahk  1 °1®S  9s[F(Sp,tq)]. 

Therefore  the  estimated  elevations  at  the  observation  grid  points  are  expressible  in 
terms  of  the  elevation  estimates  F(sp,tq)  at  the  node  points  by 

[F(x,  ,yj)]  — Ax  ox  A2  o2  (Bx  ox  B2  o2  [F(sp,tq)]) 

= (A:  Bx  ) ox  (A2  B2  ) Og  [ F (sp  ,t,)]  . 

Moreover  the  array  [F(sp  ,tq)]  is  expressible  in  terms  of  the  array  U of  observed  elevations 
by  the  equations 

t F (s  p , tq ) ] - Bq  o1  B2  o2  [ahk  1 

■t  L 

- B1  ox  B2  o2  A;  ox  A2  o2  U 
= (Bx  Aq  ) o;  (Ba  A2)  o2  U 
= (Ax  B~Y  ox  (As  B2'V  o2  U. 

This  shows  that  the  array  [F(sp,tq)J  may  be  used  as  an  array  of  coefficients  of  a functional 
model  for  the  terrain,  and  that  this  model  is  obtained  from  the  original  functional  model  by 
means  of  simple  array  operations.  The  new  functional  model  and  the  computational 
techniques  described  below  are  together  called  direct  prediction. 

As  explained  in  subsection  3.3,  the  grid  of  node  points  was  selected  from  the  grid 
of  observation  points  by  deleting  every  other  row  and  every  other  column  of  observation 
points.  Let  the  epipolar  coordinate  system  be  normalized  so  that  the  node  points  have 
consecutive  integer  coordinates.  This  requires  that  the  observation  points  be  equally 
spaced  and  form  a rectangular  grid  in  the  epipolar  plane.  Then  the  observation  points 
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themselves  include  the  node  points  plus  sets  of  points  of  the  form  (n+jy  , m),(n,m+^)  and 

(n+5,  nrl-g  ),  where  m and  n are  integers.  Once  the  values  of  N ,M  ,H  ,K , and  the  common 

_i 

values  c and  d of  the  ch  and  dk  respectively  are  assigned,  the  matrices  ,A2  and 

_i 

B2  may  be  computed  once  and  for  all. 

Since  the  values  of  F at  the  node  points  are  already  explicitly  known,  being 
the  coefficients  of  the  new  functional  model,  it  is  only  necessary,  for  the  recovery  of 
the  elevations  at  the  observation  grid  points,  to  do  computations  for  the  points  in  the 
observation  grid  of  the  form  (n,mf^),(n+^,m),(n+5,m+-^).  Accordingly  let  E,  and  E2 

_1  _i 

be  the  matrices  [fh  (rtf- f)]  Bj  and  [gk  (rrtf-|)l  B2  respectively. 

The  matrices  Ej  and  E2  have  the  following  important  property:  they  may  be 

approximated  very  closely  by  matrices  E(  and  E2  which  are  banded  and  circulant,  as 

defined  below,  except  for  the  first  and  last  few  rows.  These  properties  are  inherited 
by  and  E2  from  the  matrices 

[fh(n  + 5)]  = [exp  (c  (n  +5  - h)3)]  , 

[gk  fo+|)]  = [exp  (d  (m  + 5 -k)2)], 

Bl  = [exp  (c  (n  -h)2)], 

and 

B2  = [exp  (d  (m- k)2)]  , 

where  sh  and  tk  have  been  set  equal  to  h and  k respectively.  To  say  that  a matrix  is 
circulant  is  to  say  that  each  row,  except  for  the  top,  is  obtained  from  the  row  above  it 
by  a one  column  right  translation.  Because  their  elements  are  functions  of  the  difference 
of  the  row  and  column  indices,  the  matrices  listed  above  have  this  property. 

Because  of  their  special  structure,  E^  and  E2  are  each  essentially  determined 
by  a single  row,  and  the  only  nonzero  elements  of  the  row  are  the  ones  within  the  band. 
Finally,  symmetry  properties  of  the  matrices  result  in  a row  having  the  form 
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O...Oeb  eb_1  ...eseieie3  ...  eb_1  eb  0...0, 

where  2b  is  the  bandwidth.  These  rows  are  called  k-vectors . In  a matrix  or  array 
multiplication  involving  a matrix  E[  or  Eg , a k-vector  and  a column  would  combine 
into  an  expression 

ei(zi,J+  z,  + i/3)+  eS(2,-i,j  + 2,+2/3)+  ... 

..«+  eb(zi-b+i,j  + zi+b,3)/ 

requiring  b scalar  multiplications. 

The  sequence  |e1|,|  es  | , . . .,|  eb|  is  a monotonically  decreasing  sequence  of 
positive  numbers.  The  value  of  b is  chosen  in  the  following  manner.  Let  J be  the  flying 
altitude,  R the  base-to-height  ratio,  f the  focal  length  of  the  camera.  Also  let  the 
maximum  height  difference  from  the  mean  value  of  2b  consecutive  observation  points 
be  less  than  p%  of  the  flying  altitude  J . Finally,  let  q be  the  mean  error  of  image 
coordinates.  Then  b must  be  such  that 

Ieb+i  | (P/100) J <(J/f)R_1q. 

Typical  values  for  these  quantities  are  f = 15cm.,  R = .625  (with  image  area  equal  to 
23x23  sq.  cm.),  p = 5%,  q = 10/im.  The  quantity  J cancels  in  the  above  inequality, 
which  then  becomes  |eb+1  | < 0.02.  To  be  on  the  safe  side  and  to  allow  for  round-off 
errors,  0.02  should  perhaps  be  replaced  by  0.005. 

The  following  table  expresses  the  values  of  the  k-vector  sequence  e:  ,...,eb 
as  a function  of  the  values  of  exp(cx3)  at  x = 1 . Computations  were  performed  by  means 
of  a banded  array  Cholesky  algorithm. 
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After  the  first  few  terms,  as  for  the  other  values  of  exp(c),  the  sequence  is  very 
nearly  a geometric  progression,  with  factor  exp(c) . This  means  that  in  a matrix  or 
array  multiplication  containing  a matrix  EJ  or  Eg  with  the  above  k-vector,  the 
k-vector  and  a column  would  combine  into  dn  expression,  when  exp(c)  = 2 , 

.624(z,  + z1+1)  - . 1783(z,_1  + z,  +2) 

+ .0805(z1_s  + z!+  3)  + .004876  [z(_~ 

+z1+7  - 2[zj_5  + z1+6  + 2[-(z,_4  + z,  + 5) 

+ 2(z,_3  + z1+4)]]  - 2 [z,_7  + 

_1  .1 

zi+s+2  L -(z,_e  + z,  ^ ) + 2 (z,_g  + z1+l0)j]]. 

Here  the  second  index  j of  z has  been  supressed  to  simplify  the  appearance  of  the  formula. 

The  k-vector  term  of  smallest  magnitude  included  in  the  above  sum  is 
e10  = -.000609. 

The  evaluation  of  the  formula  requires  only  four  scalar  multiplications.  If  the 
term  e3  = .0805  is  included  in  the  approximate  geometric  progression,  only  three  scalar 
multiplications  are  needed. 

In  the  densification  , i .e. , the  recovery  of  the  elevations  at  the  observation 
grid  points  from  node  point  values  (n,m),  the  total  number  of  scalar  multiplications 
required  for  the  observation  grid  points  of  the  form  (n+§,m)  and  (n,m-t-§)  are  3(HK)  + 

3(HK)  = 6HK.  Using  one  set  of  the  elevations  so  computed,  the  elevations  at  the 
observation  grid  points  of  the  form  (n+g  ,m+g)  are  computed  with  3HK  additional  scalar 
multiplications.  The  densification  of  HK  node  point  values  by  3HK  elevations  requires 
9HK  scalar  multiplications  or  3 multiplications  per  point.  Equally  rigorous  conventional 
computational  methods  for  this  densification  require  on  the  order  of  (HK)'’  scalar  multi- 
plications . The  accuracy  required  for  the  operation  needs  in  practice  to  be  to  only 
three  or  four  decimal  places. 
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Let  bx  be  the  bandwidth  in  the  x-direction,  b2  in  the  y-direction.  The 
one-sweep  principle  of  saving  extra  read-in-read-out  operations  requires  (2bT+3)K 
words  of  core  space  for  saving  simultaneously  2bj  grid  lines  of  node  point  values  and 
3K  values  of  predictions  F(n+g,m),  F(n,m+g)  and  F(n+g  ,rrH-|).  Then  read-in  one  row 
(or  column  if  the  sweep  direction  is  changed)  of  the  grid  elevations  only  once,  and  the 
intermediate  interpolations  need  not  be  output  and  again  input  if  one  proceeds  as 
fol  lows . 


1 .  Read  in  the  elevations  of  the  (n+b^th  grid  row 


2.  Compute 

F(n+i,m) 

3.  Using  the  resulting  values  F(rr+^,m)  compute 

F(n-Ht,m+^),  m = bg  ,b2  + 1 , . . . ,K-b2 

4.  To  get  the  values  F(n,m+§),  we  use  any  of  the  2bT  rows  of  elevations 
preferably  a 

F (n+bx  , m+^ ) 

5.  Now  the  following  rows  are  no  longer  needed  in  the  computations 
and  may  be  output: 

The  (n-L)th  row  of  node  point  values 

F(n+§ , m+5 ) 

F(rrt  bx  , m+|) 

F(n+i?  ,m) 

6.  Replace  n by  n+1 . Repeat  procedure  until  n = H . 


Notice  that  the  last  interpolated  row  of  the  F(n-t-g  ,m)  and  the  F(rH-f  ,m+§)  is  for  n-N-bj  . 


The  last  interpolated  row  of  the  F(n,m+g)  is  for  n-N-b1  . 
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Because  the  node  point  elevations  and  the  interpolated  values  only  have 
5 decimals  (for  computations  only  3 or  4 decimals  are  needed  if  their  local  mean  value 
is  considered  as  a trend  function)  the  stored  2b;  K words  can  be  compressed  to  fewer 
but  longer  words  on  the  general  purpose  computers,  Gennery  [1].  If  extra  core  space 
were  available,  in  addition  to  the  necessary  2ba  K words,  one  could  apply  the  above 
algorithm  to  batches  of,  say,  p rows.  Then  the  number  of  read-out,  read-in  operations 
during  the  one-sweep  processing  is  decreased  p-fold.  The  array  processor  is  an  ideal 
development  in  this  respect  because  now  the  p rows  can  be  processed  in  parallel  as  if 
only  one  row  existed.  Especially  the  more  efficient  read-in  - read-out  operations  of 
the  parallel  processor  speeds  the  processing  compared  with  the  sequential  processor. 
Some  closer  comparisons  with  actual  processor  times  of  the  two  different  methods  was 
not  performed  during  this  study. 

The  principle  of  direct  prediction  in  the  two-dimensional  case  requires  2bl+2b3 
node  point  values  F(n,m)  only  in  the  closest  area  of  F(n+^  , nH-j),  F(n,nrt-f ),  and 
F(n+-g  ,m).  This  fact  allows  the  boundaries  of  the  area  to  assume  arbitrary  shapes  such  as 


and  the  direct  predictions  still  can  be  utilized  inside  the  narrow  zone  of  boundaries. 


4.  ARRAY  TRANSLOCATION  IN  MAP  COORDINATES  TO  PREPARE  FOR 
ARRAY  ALGEBRA  DATA  COMPACTION 

In  the  map  plane  versus  the  epipolar  plane,  the  points  at  which  elevations 
ore  observed  by  the  epipolar  scanner  do  not  form  a grid.  In  the  present  section,  a method 
called  array  translocation  is  described  for  determining  elevations  at  a set  of  points  called 
the  elevation  grid . The  elevation  grid  is  a bounded  grid  in  the  map  plane,  the  elevations 
at  whose  points  can  be  input  to  a contour  plotter.  The  set  of  elevations  at  the  points 
of  the  elevation  grid  may  be  compacted  by  the  array  algebra  method  of  Section  3,  the 
elevations  being  treated  for  this  purpose  as  an  array  of  observations  at  the  points  of  a 
bounded  grid.  Subsection  4.1  describes  the  principle  of  array  translocation.  Subsection 
4.2  contains  examples  and  algorithms  for  different  forms  of  array  translocation. 

4.1  The  Principle  of  Array  Translocation 

In  the  present  discussion,  a set  of  points  is  specified  which  in  the  map  coordinate 
system  form  a bounded  grid.  The  map  coordinate  system  is  taken  to  be  normalized  so  that 
ihegrid  is  rectangular  and  the  grid  points  have  integer  coordinates,  consecutive  grid  points 
on  a grid  line  having  consecutive  integer  coordinates.  This  bounded  grid  is  called  the 
elevation  grid.  The  problem  is  to  compute  estimates  of  elevations  at  the  points  of  the 
elevation  grid  from  a set  of  elevations  observed  by  the  epipolar  scanner  at  points  which 
do  not  form  a grid  in  the  map  coordinate  system. 

The  transformation  from  the  epipolar  plane  to  the  map  plane  is  assumed  to  be 
such  that  locally  the  departure  from  linearity  is  not  too  severe,  and  distances  are  nearly 
preserved.  Because  the  observation  points  of  the  epipolar  scanner  form  a grid  in  the 
epipolar  plane,  taken  here  to  be  rectangular,  these  observation  points  have  locally  a 
nearly  gridded  arrangement  in  the  map  plane. 
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The  method  of  estimating  elevations  at  the  points  of  the  elevation  grid  consists 
in  locally  interpolating  or  fitting  low  degree  polynomials  through  the  elevations  at  the 
observation  points,  and  then  evaluating  these  polynomials  at  the  points  of  the  elevation 
grid;  a given  polynomial  is  evaluated  at  those  points  of  the  elevation  grid  which  are 
within  the  local  region  over  which  that  polynomial  has  been  fitted. 

Any  estimation  technique  of  the  above  sort  is  called  translocation . Array 
translocation  consists  of  performing  the  estimations  in  two  steps.  The  first  step  is  an 
estimation  of  elevations  at  a set  of  intermediate  points,  with  integer  X-coordinates, 
by  translocation  involving  polynomials  of  the  X-variable  only.  The  second  step  is 
an  estimation  of  the  elevations  at  the  points  of  the  elevation  grid,  from  the  elevations 
at  the  intermediate  points,  by  translocation  involving  polynomials  of  the  Y-variable 
only. 

The  first  step  of  the  array  translocation  is  performed  as  follows.  A positive 
integer  N is  specified  in  the  interval  from  1 to  5.  The  number  N-l  will  be  the  degree 
of  the  polynomials  in  the  one  variable  X used  for  the  local  elevation  estimate  calculations. 
A second  positive  integer  M is  specified  to  be  the  number  of  consecutive  points  in  an 
epipolar  scan  line  to  which  a polynomial  of  degree  N-l  will  be  fitted.  Therefore  in 
practice,  N5M.  It  is  assumed  for  this  discussion  that  the  epipolar  scanner  scans  along 
lines  parallel  to  the  X-axis  in  the  epipolar  plane,  and  hence  in  the  map  plane  along 
curves  roughly  parallel  to  the  map  coordinate  system  X-axis. 

Let  Pk  = (xk  ,yk),  1 -k  ^M,  be  consecutive  points  on  an  epipolar  scan  lines, 
with  the  displayed  coordinates  being  relative  to  the  map  coordinate  system.  Let  zk  , 
l^kSM,  be  the  elevation  at  Pk  . Two  polynomials  F and  G of  degree  N-l  are 
required.  The  polynomial  F is  fitted  to  the  Y-coordinates  yk  at  the  X-coordinates  xk  . 

The  polynomial  G is  fitted  to  the  elevations  zv  at  the  X-coordinates  xy  . The  differences 
*<+  i ~ xk  , 1 — ic  SM-1 , are  by  assumption  all  nearly  equal  to  1 . For  the  purpose  of 
determining  the  coefficients  of  F and  G,  these  differences  will  be  taken  to  be  equal  to  1 . 
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It  is  known  from  the  adjustment  calculus  that  this  approximation  will  not  introduce 
any  serious  errors  in  the  estimation.  Let  H be  the  largest  integer  less  than  or  equal  to 
M/2.  Then  H+h,  with  h assuming  integer  values  in  the  interval  -H  < h S(M+l)/2, 
will  be  the  numbers  replacing  the  X-coordinates  x1,.../xH  for  the  purpose  of  fitting 
the  polynomials.  For  the  purpose  of  the  fittings,  all  X-coordinates  will  be  translated 
by  -H  , so  that  X-coordinates  of  the  points  to  which  the  polynomials  are  fitted  are 
taken  to  be  the  integers  h,  -H5  hS  (M+ 1)/2 . Set 

F(x)  = L a,  x*  , 

1 = 0 

and 

G(x)  = Y b.  x'  . 

1=  o 

Let  C be  the  (MxN)  matrix 

C « Chi'1]. 

T T 

where  hk  = -H  +k  , 1 § M , and  1 St  5 N . Let  Y = [y1  . . .yM  ] and  Z = [z1  . . .z.,  ] . 

Also  let  A = [q,  . . .aN  ] T and  B = [t^  ...bN_x]T.  Then  the  values  of  A and  B for 

which  the  polynomials  F and  G give  the  best  fit,  in  the  sense  of  least  squares,  to  the 
respective  sets  of  pairs  of  numbers  (hk  ,yk  ) and  hk  ,zk  ),  lSk  SM,  are  given  by 

A = Yy 

and 

B = ClZ. 

> 

Here,  as  in  Section  2,  C = (CT  C)  CT  . 

I 

It  is  clear  that  for  fixed  N and  M,  that  the  matrices  C and  C are  independent 
of  the  particu  lar  sequence  of  points  Px  , . . . , P„  , and  therefore  in  an  algorithm  for 
computing  the  column  vectors  A and  B at  the  various  intervals  along  the  scan  lines, 

C may  be  specified  once  and  for  all. 
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The  polynomials  F and  G,  with  their  coefficients  determined  as  above,  are 
used  to  obtain  elevation  estimates  at  points  with  integer  X-coordinate  as  follows. 

Let  m be  the  integer  closest  to  the  midpoint  of  the  interval  from  xx  to  . Then 
G(m-H)  is  the  estimated  elevation  at  the  point  (m,F(m-H)). 

The  first  step  of  the  array  translocation  is  complete  when  all  consecutive 
sets  of  M observation  points  on  the  epipolar  scan  lines  have  been  processed  in  the 
above  fashion  to  produce  a new  collection  of  points  all  having  integer  X-coordinates 
and  estimated  eievations. 

The  second  step  is  a repetition  of  the  first  step,  but  now  along  lines  parallel 
to  the  map  coordinate  system  Y-axis  and  intercepting  the  map  coordinate  system  X-axis 
in  points  with  integer  X-coordinates.  The  points  to  which  the  polynomials  are  fitted 
are  the  ones  produced  in  the  first  step;  the  points  at  which  the  elevations  are  estimated 
are  the  ones  with  both  coordinates  integers.  These  are  the  points  of  the  elevation  grid. 

This  completes  the  general  description  of  array  translocation.  Given  below  are 
some  examples  of  the  basic  translocation  fitting  and  estimation  procedures  for  various 
values  of  N and  M,  accompanied  by  efficient  algorithms  for  executing  the  computations. 

4.2  Examples  and  Algorithms  for  Array  Translocation 
a)  N = 1,  M = 1 

For  this  case  the  elevation  estimate  at  a point  Q with  integer  X-coordinate 
m is  the  value  of  the  elevation  at  the  epipolar  scanner  observation  point  nearest  to  Q. 

The  Y-coordinate  of  Q is  the  Y-coordinate  of  this  observation  point. 
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This  is  ordinary  linear  interpolation. 


c)  N = 2,  M = 3 


Here 


C ' = 


o h 


and  Y and  Z have  size  (2x1). 

Then  if  m1  = m-H  , 

F(m')  = [1  m']A 

= [1  m']cV 


An  algorithm  for  evaluation  F at  m‘  is  as  follows: 

di  = (y»  + y2  + y3)/3 
ds  = (y3  - yi) 2 

F(m')  = dj  + d2  m1 . 

This  requires  only  two  multiplications  by  the  constants  3 and  2 , and  one  multiplication 

d m'.  The  function  G is  evaluated  at  m*  by  a similar  algorithm  with  Y replaced  by  Z. 


d)  N = 2,  M = 4 

For  this  case, 


C 


l 


6 


2 1 0 
-1  13 


and  Y and  Z have  size  (4x  1). 


An  algorithm  to  evaluate 

F(m')  = [1  m']  C^Y 


is  as  follows: 

di  = yi  + ys 

d2  = dl  +dx 

d3  = Yi  + y3 

d4  = (dL  +ds)/6 

ds  = (y4  - y s) 2 

4s  = d4  + ds 


F(m')  = d5  + 4 m'. 


Even  though  there  are  now  four  observations,  evaluation  of  each  of  the  polynomials  still 

_i  -1 

only  require  two  multiplications  by  the  constants  2 and  6 , and  the  multiplication 

<4  m1,  which  is  the  same  number  as  is  required  in  case  (c)  where  M = 3.  The  fourth 
observation  may  cause  some  oversmoothing;  however,  including  it  improves  the  chances 

•4 

for  the  detection  of  gross  errors.  This  may  be  done  by  forming  L z®  - 4G(m'),  and  if 
it  is  unexpectedly  large,  the  residuals  z,  - G(x, ) may  then  be  tested  to  determine  the 
point  P{  at  which  the  gross  error  occurs. 


e)  N = 3,  M = 5 
Here 


-3/35 

12/35 

17/35 

12/35 

-3/35 

-1/5 

- 1/10 

0 

1/10 

1/5 

1/7 

- 1/14 

- 1/7 

1/14 

1/7 

and  Y and  Z have  size  (5x  1). 
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Instead  of  estimating  the  elevation  at  a single  point,  simultaneous 
estimation  at  three  points  will  be  considered  in  this  example.  These  three  points  are 
(m,/F(m1))  and  also  (m*-1 , F(m'- 1))  and  (m’+  1 , F(m'+ 1)) . 

The  algorithm  for  evaluation  of  G at  m'-l,  m1  and  m'+l  is  as  follows: 

dj  = zs  + z4 

ds  = dx  + dx 

d3  = da  + ds 

d4  = Zl  + *5 
d5  = d!  - Z3 
4 = d4  + 4 

d7  = (3/35)(-d4  + d3  + (17/3)zs) 

d8  = (l/10)(z4  - z2  + zs  - zx  + Zg  - zT) 
ds  = (1/14)(4  -dx). 

So  far,  only  four  multiplications,  each  with  one  factor  fixed,  are  required.  Then 

e,  = <4  m' 

es  = 4 + ei  + e! 

F(m')  = d7  + m'(dq  + e:) 

F(m'+1)  = F(m')  + dg  + es 
F(m'-l)  = F(m')  + d,  - es  . 

Two  more  scalar  multiplications  are  required  above. 
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This  case  has  the  advantage  of  not  oversmoothing,  as  well  as  producing 
three  elevation  estimations  at  each  iteration. 

If  a sparser  elevation  grid  is  wanted,  e.g.,  only  at  points  with  even 
integer  map  coordinates,  then  the  last  part  of  the  above  algorithm  may  be  easily 
modified  for  this  as  follows: 


e! 

d9m' 

e2 

dq  + ex  + ex 

e3 

d9  + dg 

e4 

e2  + e2 

% = 

e3  e3 

II 

LI- 

d7  + m‘(dB  + ex) 

F(m'+2)  = 

F(m')  + e4  + e$ 

F(m'-2)  = 

F(m')  + e4  - e5  . 
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5. 


SUMMARY  AND  CONCLUSIONS 


In  the  previous  sections  the  principles  of  array  algebra,  as  it  applies  to 
digital  terrain  representation,  were  put  forth.  Section  2 detailed  the  mathematical 
formalization  for  array  regression  in  an  N-dimensional  Euclidean  space.  Section  3 
put  forth  the  algorithms  for  data  compaction  in  an  epipolar  coordinate  frame  while 
Section  4 detailed  the  methods  of  array  translocation  as  it  would  be  applied  to 
generating  rectangu  larly  sea  led  terrain  elevation. 

During  this  study  it  was  found  that  efficient  direct  array  prediction  algorithms 
could  be  derived  and  that  the  explicit  formulation  of  the  regression  coefficients  could 
be  derived  and  that  the  explicit  formulation  of  the  regression  coefficients  could  be 
avoided.  Using  these  algorithms  the  array  requirement  for  rectangularly  spaced  data 
was  removed;  thereby  allowing  boundaries  of  the  area  to  assume  arbitrary  shapes.  In 
delineating  these  methods  several  examples  of  direct  prediction  coefficients,  or  "k- 
vectors",  were  derived.  In  typical  cases,  the  knowledge  of  ten  (10)  such  coefficients 
replaces  the  data  handling  as  has  been  conventionally  required  for  explicit  formulation 
of  observation  and  normal  equations,  solution  of  normal  equations  and  predictions  using 
the  solution . 


The  algorithms  developed  for  translocation  demonstrate  that  between  2 and  10 
multiplications  are  required  for  each  point  manipulated. 

In  addition  to  analyzing  the  mathematical  equations  required  for  terrain 
compaction,  the  computational  requirements  were  analyzed  for  both  sequential  and 
parallel  processors.  Based  upon  this  analysis  it  appears  that  the  array  processor  will 
provide  a significant  computational  advantage. 

As  was  previously  stated  two  alternate  methods;  namely,  compaction  and 
translocation,  were  reviewed  and  mathematically  described  herein.  Based  upon  our 
analysis  it  appears  that  these  processes  could  have  a significant  impact  on  the 
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computational  efficiency  in  the  field  of  automated  cartography.  For  this  reason  we 
recommend  that  these  algorithms  be  implemented  at  the  ETL  facility  and  a series  of 
comparative  tests  conducted  using  both  the  CDC  6400  and  Staran  computer  systems. 

The  exploitation  of  Staran,  in  conjunction  with  array  regression,  should  form  the 
cornerstone  for  a technology  for  efficient  manipulation  of  massive  volumes  of  digital 
terrain  data.  Increases  in  overall  throughput  of  orders  of  magnitudes  are  totally  within 
current  state-of-the-art. 
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APPENDIX  A 

REVIEW  OF  THE  R-MATRIX  MULTIPLICATION  PRINCIPLE  OF  ARRAY  ALGEBRA 
General  Background  and  Definitions 

This  expose  is  intended  to  explain  some  technical  details  pertinent 
to  the  R-matrix  multiplication  which  forms  the  backbone  of  the  so-called 
Array  Algebra,  treated  in  [Rauhala,  1976]  and  in  many  other  publications  by 
the  same  author.  Considering  the  limited  scope  of  this  paper,  the  proofs  of 
equivalence  between  certain  least  squares  problems  solved  in  a conventional 
way  and  the  extremely  efficient  solutions  of  the  same  problems  via  array 
techniques  will  be  omitted.  We  shall  be  satisfied  by  observing  that  the  con- 
ventional (monol inear)  least  squares  solution  may  be  generalized  to  multilinear 
cases  of  array  algebra  provided  the  design  matrix  A in  the  monol inear  formula- 
tion can  be  expressed  as  a tensor  product  of  several  matrices, 

ft  = ft,  ® ftz  ® . . • ® , ( 1 ) 

(m-n)  ( «,■«,)  ('VnJ 

M * rrt,  » mz  * . . • « mi  , > / , 

N = n,  * n2  - .. . « r?i  . J 

Here  the  symbol  ® represents  a tensor  product  and  we  are  in  the  presence  of 

an  i-dimensional  case.  The  tensor  product  itself  is  defined  in  terms  of  two 
matrices  as  follows: 

B ® C * IB  (C)„j  , 

(r « i)  (Pr*?4) 

i.c. , each  "entry"  of  the  resulting  matrix  ic  the  entry  of  the  second  matrix 
multipl ied  by  the  whole  first  matrix.  From  this  definition  follows  e.g.  the 
computational  rule 

3®  C ® J)  s £ =(^®C)®J)=  3 ® (C  ® J)  ) . 

(P’<j)  f r.S)  (<.u)  (prt.^iu)  (pr  . (f*u;  (P’j)  ( rt’iU ) 

We  only  mention  that  a formulation  such  as  (1)  is  possible  if  the  observations 
form  an  i-dimensional  grid;  throughout  this  paper,  the  pertinent  weight 
matrices  will  be  considered  as  unit  matrices  for  simplicity. 


The  fami  1 iar  observation  equations  in  the  mono-linear  case  read  as 
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f 

L 
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follows: 

V = P X - l , 

( M ■ I ) (M  * *0(N  ■ O ( M ■ I ) 

where  v is  the  residual  vector,  A is  the  design  matrix,  x is  the  vector  of 
unknown  parameters,  and  t is  the  vector  of  observations  (if  the  problem 
were  not  linear  at  the  outset  it  would  be  linearized  and  a constant  vector 
would  be  added  to  the  right-hand  side).  The  least  squares  solution  may  be 
written  as 

(x).  = I {Al)jr(t)r  , j-  1.2,  ....  N , (2a) 

where  (A*  )^r  represents  the  entries  of  the  matrix  A*  and  (x)^  , ( l)r 
represent  the  entries  of  the  vectors  x ,1  respectively; 

)"  At 

( N«  M)  (N • H)  ( N * M ) 

is  called  the  1-inverse  of  the  matrix  A under  the  usual  assumption  of  the 
nonsingular  matrix  of  normal  equations,  ATA.  Written  in  the  matrix  form, 
equation  (2a)  becomes 

x = ft1  l ■ (21) 

( N - I ) (N « M ) (M  * I ) 

We  shall  now  confine  our  attention  to  an  array  (multilinear)  form- 
ulation associated  with  (1).  We  notice  that  t will  be  written  as  an 
(m( « m2- ...«  m^)  array  L which  stems  from  the  gridded  observational  structure. 

In  a similar  fashion,  x will  be  represented  by  an  (n,«  n2*...  >04.)  array 
X.  Furthermore,  the  matrix  A*  can  be  shown  to  be  the  tensor  product 

A1  = Alt  ® A[  ® ...  ® A[  . (-2c) 

(N«M)  (n2«m2)  l°i  *mi) 
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From  various  references  on  the  subject  of  array  algebra  we  observe  that 


the  solution  equivalent  to  (2a),  (2b)  is  the  following: 

(XL..,  -lmrL  fftfULU 

J,  = I , Z } . • • ; n, 

jz  ~ ^ ^ > n2 


ji  ~ I > ^ ‘ > ni 

In  a compact  form,  this  may  be  written  as 


X = 

( •V'V  * nJ 


(A,*)'  (AH 

(",«•",)  (»»*"'*) 


t \2 


(fl‘r  l . 

(nt-m.)  . nrt) 


(31) 


An  explanation  of  this  notation  convention  is  in  order.  The  symbols  X and  L 
represent  i-dimensional  arrays,  while  Aj  , j=l,  2,...,  i,  are  matrices  (two 
dimensional  arrays).  The  subscript  simply  identifies  the  matrix,  i.e.,  it 
has  a role  of  a "name".  Thus  instead  of  A,  , A2>  A3  , etc. , we  could  have 
written  A,  B,  C,  etc.;  since  a matrix  A^  is  assumed  to  have  the  dimensions 
(mj»nj),  its  1-inverse  matrix  A*  has  the  dimensions  (n^«m^)-  On  the  other 
hand,  the  superscript  "k"  of  a matrix  indicates  that  its  second  index  corresponds 
to  the  k-th  index  of  L.  For  instance,  the  matrix  A*  has  been  assigned 
the  superscript  "i";  this  means  that  the  second  index  of  the  elements  in 
A^  must  be  the  same  as  the  i-th  index  of  the  elements  in  L.  Equation  (3a) 
demonstrates  this  rule  and  shows  that  the  index  in  question  is  r . Consequently, 
the  order  of  the  second  dimension  of  Af  in  (3b)  is  the  same  as  that  of  the 
i-th  dimension  of  L,  namely  mt  . In  analogy  to  the  conformabi 1 ity  requirement 
(the  number  of  columns  in  a first  matrix  equals  the  number  of  rows  in  a 
second  matrix)  associated  with  the  conventional  matrix  multiplications,  we 
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may  now  speak  of  conformable  arrays  associated  with  the  R-matrix  multiplications. 
The  right-hand  side  of  (3b)  is  a typical  example  of  i successive  R-matrix 
multiplications  in  a convenient,  "short-hand"  notation.  A "long-hand" 
equivalent  of  (3b)  is  of  course  equation  (3a)  in  which  j,  , are 

the  "free"  indeces.  It  is  worth  mentioning  that  the  order  in  which  the 
matrices  are  arranged  for  the  R-matrix  multiplications  is  inconsequential 
provided  the  superscripts  are  not  interchanged;  equation  (3b)  could  thus 
be  written  e.g.  as 

x = {Atf  [ft!)' ...  (W  L . 

This  is  apparent  from  (3a)  in  which  the  summations  may  be  performed  in  an 
arbitrary  order. 


We  shall  next  illustrate  a few  special  cases  of  R-matrix  multiplications 
If  the  space  is  one-dimensional  (i=l),  we  have  essentially 


X = (/),')'  L * fl/L 

(VU  (W,«r »,)  ("V*) 


which  corresponds  to  (2b).  The  two  dimensional  case  may  be  portrayed  as 


X - (A/)1  (flin  -=  /),*  L R" 


where  A^T=  (a|  )t  . That  it  is  so  becomes  apparent  from  the  "long-hand" 
form,  namely 


(x). . - £(*,«).,  £rfl/u Drr  -=  zr  «l(l)(  Kr 

V 'j,)L  r,.i  1 J>r>  rr-,  2 J*r*  r.ri  r,x.  rA.|  }<r,  r.r, 

where  (A*  ).  f has  been  written  as  (A$ T)r^  . We  can  thus  assert  that  the 
superscript  "1"  corresponds  to  the  left-hand  side  multiplication  and  the 


superscript  "2"  corresponds  to  the  transposition  and  the  right-hand  side 
multiplication.  Unfortunately,  there  is  no  equivalent  for  the  superscript 
"3",  etc.,  so  that  somewhat  absurd  names  such  as  "back-hand"  side  or 
"beyond-hand"  side  multiplication  would  have  to  be  invented  if  we  wanted 
to  illustrate,  at  least  intuitively,  some  operations  in  terms  of  the  conven- 
tional matrices.  One  might  thus  be  tempted  to  visualize  the  R-matrix  multi- 


plications in  three  or  i dimensions  as  follows: 

X - CA1)'  M/r  L = 

(V',z*nj)  (v^,)  ('*>, 


ft!  L A[ 


4 1 


X = (Alt)'(Ait)1(A'),...(A!)iL  * a!  L H"  . 

However,  it  should  be  sufficient  to  work  solely  in  terms  of  the  original 
definitions  (3a)  and  (3b).  We  shall  next  present  some  practical  consider- 
ations in  handling  equation  (3a)  and  shall  determine  the  number  of  significant 
computer  operations  (scalar  multiplications)  involved. 

The  form  of  equation  (3a)  suggests  that  each  summation  sign 
represents  a "do-loop"  in  regular  FORTRAN  programming.  The  last  do-loop 
may  then  be  imagined  to  result  in  an  array  Z whose  elements  are 

(ZL...r  , - £ . (U • <3t) 

In  computing  the  (one)  element  of  1 depicted  in  (3c),  m^  multiplications 
must  be  performed  with  regard  to  the  index  r which  is  then  replaced  by  the 
free  index  . To  fill-in  the  elements  of  Z for  one  value  of  , the 
summation  in  (3c)  which  was  seen  to  require  n\  multiplications  must  be 
carried  out  m,  * m^.-.-m^.,  times  (m,  times  due  to  the  index  r(  , m^  times 
due  to  r2  , etc.).  Finally,  such  operations  must  be  repreated  for  each  , 
i.e.  times,  in  order  to  complete  the  array  Z.  This  step  has  accordingly 
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invol ved 


m(  < ( * ..."  m.ml  ) * ni 


multiplications.  In  analogy  to  (3a)  and  (3b)  we  may  now  write 

'%  v.  ■ !<**.  !>•%« 

x = (/),')■  K)1...  (fl',r  z 

K«"1,)  frti-rrJi)  (ni-l  * mM  ) («,«  "Hi ) 

A similar  argument  to  the  above  may  be  repeated  with  respect  to  the  last 

do-loop  in  (3d)  giving  rise  to  the  elements  (Y)  of  , ^ . 

which  now  involves 


"V,  - ( m,  « « . . . - "V*  « ni  ) « n^, 

multiplications.  The  number  of  multiplications  that  occur  in  the  remainder 

of  the  i steps  may  be  similarly  expressed  as 

m^2  * ) " ni-2 


m,  * ( * rij  * . . . * n i ) * n,  . 

It  may  be  appreciated  that  although  the  result  is  invariant  with  respect  to 
the  order  in  which  the  array  multiplications  are  carried  out,  the  total 
number  of  scalar  multiplications  is  not. 

In  order  to  facilitate  the  upcoming  comparisons  in  terms  of  scalar 
multiplications  (and  thus  in  terms  of  the  computer  time),  the  following 


simplifications  are  made: 


i,  = m2  = . . . = mt  * m , 


n,  » n2  = ...  = ni  - n 


[Aa) 

(4*) 


-6  - 


• i 


we  have 


m = k n (5a) 

where  necessarily 

k } \ . (St) 

The  number  of  scalar  multiplications  associated  with  the  R-matrix  multiplications 

as  described  in  the  preceding  paragraph  is  now 

m * ( m*'"1 ) « n = n"  k* 
m « (rn  2 n)*  n - k 
m * (m  3 n1)' n = n' 1 kl 
# 

m * ( ri1  ' ) * n = ^ k . 

The  total  is 


l n . . . if  k = / > 

(7a) 

n;*’k  JLlL  ...if  k>l  . 
k - 1 

(U) 

The  formation  of  each  of  A:  from  Aj  requires  approximately  only 

!•*•■»»>  r 

(2k  +l)n3  scalar  multiplications  (mn2  = kn3  to  form  Aj  Aj , approximately  n3  to 
invert  this  positive-definite  matrix,  and  kn3  to  post-mul tiply  this  new  matrix 
by  Aj).  The  total  in  an  i-dimensional  problem  is  thus  i (2k  +1)  n3  . If  i>2, 
this  number  is  at  least  one  order  of  magnitude  smaller  than  the  number  of 
scalar  multiplications  associated  with  the  subsequent  R-matrix  multiplications 
as  evidenced  by  (7) . Under  such  circumstances,  the  computer  time  required  to 

l 

form  all  of  Aj  may  be  considered  negligible. 
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Efficiency  of  the  Array  Technique 


In  cases  where  the  design  matrix  has  the  structure  exhibited  in 
(1),  the  conventional  procedure  of  forming  the  normal  equations  and  solving 
for  x as  in  (2b) , i .e. 


x = flT  l 

(N«l)  f /V  * N ) (N‘M)  (M‘l) 

would  be  extremely  uneconomical.  In  particular,  the  inversion  of  A1" A would 
necessitate  approximately  N3  scalar  multiplications.  This  inversion  would 
consume  by  far  the  greatest  amount  of  computer  time.  The  formation  of  AT£  , 
by  comparison,  would  require  only  NM  scalar  multiplications  (although  in 
general  M>N,  M is  of  the  same  order  of  magnitude  as  N)  and  the  pre-multi- 
plication of  this  vector  by  the  above  inverse  matrix  would  necessitate  only 
Na  ?'J-'tional  scalar  multiplications.  We  have  not  mentioned  the  original 
construction  of  the  matrix  A1"  in  which  slightly  more  than  NM  multiplications 
would  be  performed.  In  any  event,  the  number  of  significant  scalar  multi- 
plications determining  the  computer  run  time  in  this  "brute  force"  approach 
would  be  N3  and  the  operation  directly  responsible  for  it  would  be  the  matrix 
inversion  of  ATA.  Before  proceeding  any  further  we  emphasize  that  the 
monolinear  solution  could  be  carried  out  much  more  efficiently  through  the 
use  of  the  property  (2c).  If  we  set  for  simplicity 
n,  = n2  = ...  = - n , 

the  "brute  force"  approach  would  involve  one  inversion  of  an  (nl*  n‘ ) matrix 
while  a "ref ined"monol inear  approach  will  be  seen  to  necessitate  i inversions 
of  small  (n«n)  matrices.  Thus  if  n = 10,  we  would  have  e.g. 

i = 2,  N = nl  = 100  parameters;  ( 100  * 100)  inverse  versus  two  (10  * 10) 

i = 3,  N = n1-  = 1000  parameters;  (1000  * 1000)  inverse  versus  three  (10  * 10) 
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The  inversions  of  (n*n)  matrices  account  for  a negligible  amount  of  the 
total  computer  time.  On  the  other  hand,  the  inversion  of  an  (N  * N)  matrix 
in  the  "brute  force"  approach  is  in  general  totally  prohibitive;  this  would 
not  be  substantially  altered  even  if  an  inverse  multiplication  of  the  vector 
At*  by  the  matrix  ATA  were  performed  without  explicitly  inverting  ATA.  The 
total  number  of  scalar  multiplications,  although  reduced,  would  still  be 
proportional  to  N3.  The  associated  computer  storage  problems,  not  treated 
in  this  paper,  would  be  equally  formidable. 

For  the  reasons  exposed  above,  only  the  "ref ined"monol inear  approach 
will  be  now  considered  and  later  compared  with  the  array  solution  for  efficiency. 
When  the  tensor  product  property  (2c)  is  utilized, from  (2b)  and  (2c)  we  write 

X = ( f\f  ® f\[  (s>  ...®  f\[  ) l . (8) 

The  formation  of  A(?  ® a|®.  . . ® a[  requires 

f * ...  + 0(rr)(  * n* rttj * . . . « nxmx  (9) 

scalar  multiplication,  all  of  A j having  been  computed  beforehand.  Due  to 
(l'),  the  last  term  in  (9)  is  MN  and  the  next  largest  term  is  smaller  by  at 
least  two  orders  of  magnitude.  Since  the  number  of  scalar  multiplications 
needed  to  form  A(t  in  (8)  is  also  MN,  their  total  number  is  then  2MN. 

Instead  of  (8),  one  might  consider  a mathematically  equivalent  hut 
formally  slightly  different  equation,  namely 

x = [(47),)"®  (/£*,)"....•  ] aTi . 

However,  in  this  formulation  substantially  more  scalar  multiplications  are 
needed.  Already  the  vector  ATf  involves  2MN  of  them  since 

$r=  fl,r®  Pi[®  .. . ® Al 

is  in  this  respect  equivalent  to  Af;  its  subsequent  post-multiplication  by  l 
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necessitates  another  MN  scalar  multiplications.  Furthermore,  the  construction 
of  the  matrix  between  the  brackets  in  the  above  equation  accounts  for  N 1 
multiplications  as  does  the  final  multiplication  by  ATf.  For  these  reasons, 

the  most  advantageous  route  pursued  in  the  refined  monol inear  formulation 
is  that  of  equation  (8). 


We  now  adopt  the  same  simplifications  as  considered  earlier  in  the 
equations  (4)  and  (5)-,  it  follows  that 

M = ~ kv  , ( 10a) 

H = ■ (lOt-) 

The  total  number  of  scalar  multiplications  associated  with  (8)  is  thus 

2 MN  = 2k' nH  . O') 

The  number  of  multiplications  used  in  forming  all  of  f\j  was  seen  following 

(7)  to  be  i (2k  +1)  n3.  This  number  is  at  least  one  order  of  magnitude 

0 

smaller  than  (11)  even  if  i = 2 and  the  formation  of  the  matrices  Aj  may 
accordingly  be  neglected  insofar  as  the  computer  time  is  concerned.  We  shall 
next  compare  the  above  most  favorable  monolinear  formulation  with  the  array 
approach  for  efficiency.  As  a matter  of  interest  we  mention  that  if  k = 1 
(i.e.  n=m,  N = M),  this  monol  inear  formulation  exhibits  the  structure  related 
to  the  Fast  Fourier  Transform. 


We  have  seen  that  both  the  refined  monolinear  formulation  and 

l 

the  array  approach  require  the  formation  of  the  Aj  matrices.  In  the  former 
case,  these  operations  are  negligible  compared  to  the  total  computer  burden 
for  any  i ^ 2 while  in  the  latter  case,  their  impact  can  be  ignored  only  if 
i>  2 (for  i = 2,  the  number  of  pertinent  scalar  multiplications  would  have 
the  same  order  of  magnitude  as  in  the  subsequent  R-matrix  multiplications). 
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Accordingly,  the  "savings  ratio"  that  will  be  determined  in  conjunction  with 
the  R-matrix  multiplications  will  represent  also  the  "total  savings  ratio" 
only  for  i > 2.  If  i = 2,  the  total  savings  ratio  of  the  array  solution 
versus  the  refined  monolinear  solution  would  be  somewhat  lower,  although 
its  order  of  magnitude  would  not  be  altered.  The  savings  ratio  (R)  that 
concerns  us  is  given  by  the  expression  in  (11)  divided  by  the  number  in  (7a) 
or  (7b).  We  thus  have 

R = (2/ i)  if  k*l  , (12a) 

R = 2ri‘-'(ki-ki")/(ki-l)  ...if  k>\  . (lit) 

For  a large  value  of  k in  (12b),  i.e.  for  many  more  observations  than  para- 
meters, it  follows  that  R-»2nv~'  . For  example,  if  i = 3 and  n = 100  resulting  in 
N = 1000000  parameters  , the  (total)  savings  ratio  would  be  R-»  20000.  However, 
a resonable  value  for  k in  practice  is  k = 2;  one  can  imagi ne  that  every  second  point 
in  each  dimension  of  the  grid  is  a "node  point"  where  the  value  of  some 
function  is  to  be  determined  from  the  least  squares  adjustment. 

We  conclude  this  section  with  a few  examples  to  illustrate  the 

savings  one  may  achieve  when  utilizing  the  array  techniques  : 

i = 2,  n = 10,  N = 100  parameters...  k = 1,  R = 10.0 

k = 2,  R = 13.3 

k = 4,  R = 16.0 

k = 8,  R = 17.8 

k large  R -*■  20.0 

i = 3,  n = 10,  N = 1000  parameters...  k = 1,  R = 67 

k = 2,  R = 114 

k - 4,  R = 152 

k = 8,  R = 175 

k large  R -*  200  . 

In  this  last  group  (i  = 3),  R has  the  role  of  the  total  savings  ratio.  The 
number  n = 10  would  of  course  be  too  small  in  practice.  A quite  realistic  example  in 
two  dimensions  might  be  the  following: 
i = 2,  n = 100  , N = 10000  parameters...  k = 2,  R = 133. 
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It  is  not  difficult  to  realize  that  if  the  "brute  force"  approach  were  used 
in  the  cases  that  concern  us  (gridded  observations),  the  savings  ratio  would 
be  several  orders  of  magnitude  higher  than  the  above  R.  In  particular,  it 
would  be  multiplied  by 

(N3  + 2MN+N2)/2MN  , 

i.e.  approximately  by 

N2/2M  = N/2kL 

if  k is  not  excessively  large.  In  practice  N (total  number  of  parameters) 
could  range  from  several  hundred  to  several  million.  This  thought  should 

■ 

merely  warn  the  "problem  inventor"  that  extremely  significant  computer 
savings  may  fail  to  materialize  if  enough  care  is  not  taken  at  the  level 
of  the  problem  design  and  formulation. 
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Future  Saving s Possibilities 

VJe  assume  that  an  ideal  parallel  processor  capable  of  performing  a 
sufficient  number  of  scalar  multiplications  simultaneously  will  become  avail- 
able in  the  future.  The  economies  achieved  in  conjunction  with  the  array 
techniques  could  then  be  greatly  magnified.  In  particular,  the  parallel 
processor  could  perform  the  multiplications  in  (3c)  involving  some  or  all 
r(  , r2,...r-  simultaneously.  In  considering  for  our  purpose  the  computer 
time  requirements  in  the  latter  case,(L)r_  _ would  act  as  a vector  with 
a simple  index.  For  the  forthcoming  comparisons  we  again  adopt  the  simplifica- 
tions (4a)  and  (4b).  At  each  of  the  i steps,  the  number  of  multiplications 
indicated  between  the  parentheses  in  the  expression  (6)  would  be  now  performed 
simultaneously  with  the  number  of  multiplications  (m)  factored  before  the 
parentheses.  The  total  number  to  be  considered  in  this  context  is  thus 

imn  = ikn1  . (|3J 

Upon  consulting  (7a)  and  (7b),  we  deduce  that  the  additional  savings  ratio 
(R)  due  to  the  parallel  processor  is 

R * n1"  ...  if  k*l  . (14a) 

R = (n"'/v  ).  (k‘-  l)/(k-  I ) ...  il  k>l  . (Ht) 

The  pertinent  savings  ratio  (R)  between  the  refined  monolinear  solution  and 
the  array  solution  in  conjunction  with  the  parallel  processor  is 

R - R R . (is) 

Upon  considering  (12)  and  (14),  for  any  k £ 1 we  have 

R - (2/i)ki"'na'1’  . 

• A • 

This  result  is  verified  if  the  number  2kc  n of  (11)  is  divided  by  ikn2  of  (13). 
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We  conclude  by  extending  some  of  the  examples  of  the  previous 

section  as  follows  (the  values  of  R are  rounded): 

i = 2,  n = 10,  N = 100  parameters...  k = 1,  R = 10  , ft  = 100 

k = 2,  R = 15  , ft  = 200 

k = 4,  R = 25  , R = 400 

i = 3,  n = 10,  N = 1000  parameters...  k = 1,  R = 100  , ft  = 6700 

k = 2,  R = 233  , R = 27000 

k = 4,  R = 700  , R = 107000 

i = 2,  n = 100,  N = 10000  parameters. . . k = 2,  ft  = 150  , ft  = 20000. 

As  a closing  remark  we  note  that  the  idea  of  a parallel  processor  suits 


perfectly  the  R-matrix  multiplications. 


Practical  Applications 


The  usefulness  of  the  R-matrix  multiplication  technique  for  a 
whole  class  of  least  squares  problems  has  been  demonstrated  at  the  level 
of  the  solution  itself.  The  rule  governing  these  multiplications  has 
been  explained  following  equation  (3b),  in  which  A*  were  the  1-inverses 
of  the  "original"  matrices  A^ . However,  essentially  the  same  rule  could 
be  applied  and  comparable  savings  could  be  realized  with  regard  to  the 
least  squares  predictions  in  which  some  "original"  matrices  would  be 
involved.  In  either  case  the  matrices  attributed  a superscript  are 
called  R-matrices.  Practical  applications  associated  with  a special 
structure  of  these  matrices  will  be  briefly  mentioned  in  the  next  paragraph. 

In  practice  the  R-matrices  often  exhibit  a banded  structure  which 
can  significantly  accelerate  the  R-matrix  multiplications.  The  banded 
structure  is  the  basis  of  an  extremely  efficient  multilinear  prediction 
technique  called  array  prediction.  According  to  [Rauhala,  1976],  this 
technique,  when  applicable,  has  the  potential  to  surpass  several  "fast 
transforms",  local  polynomials,  and  other  conventional  signal  processing 
and  modeling  techniques  in  many  sciences.  This  is  especially  true  for 
the  automation  of  the  typical  "map"  sciences  like  geodesy,  photogrammetry, 
cartography,  meteorology,  oceanography,  geology,  image  processing,  remote 
sensing,  and  other  sciences  of  the  space  age  technology  as  indicated  in 
the  above  reference. 
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